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Abstract 


The  paper  surveys  the  approximation  and  convergence  properties  of  the  h-p 
version  of  the  finite  element  method,  the  basic  theoretical  results  together 
with  the  main  ideas  of  the  proofs.  Numerical  examples  are  given. 


1. 


Introduction 


The  h-p  version  of  the  finite  element  method  is  characterized  by  the 

family  9^  =  {JR}  of  meshes  DR  covering  the  domain  Q  defining  the  elements 

J,  JA.  -  {y}  and  by  the  distribution  qCJR)  of  the  degrees  pO*)  of  the 

elements  S'  €  3R.  The  shape  functions  are  the  usual  "pull-back"  polynomials 

i.e.  mapped  polynomials  of  degree  pO")  defined  on  the  standard  element  T. 

The  degree  pi^)  can  be  different  in  different  directions  in  T.  The 

standard  element  can  be  a  triangle  or  square  in  two  dimensions  and  cube, 

wedge  or  simplex  (tetrahedron)  in  the  three  dimensions.  We  denote  by 
1 

S(3n,q)  c  H  (Q)  the  finite  element  space  and  by  N(!in,q)  its  dimension. 

The  main  problem  in  the  finite  element  method  is  to  estimate 

(1.1)  inf  |u  -  =  $(u,f,q,N). 

S  (On.q) 

N(0D  ,  q)  £  N 

DTleg: 


Assume  for  example  that  ?  is  a  (reasonable)  family  of  quasiuniform  meshes 

on  a  (reasonable)  domain  £2  c  n  =  2,3  and  that  the  degree  of  all 

elements  is  the  same  i.e.  p(3^)=p5pQ,pil,  JeJIl.  Further  assume  that 
s 

u  e  H  (n).  Then  we  have 


(1.2) 


$(U,3^,q,N)  5  llulljjs 


with  fi  =  minCp^,  s-1). 

If  we  know  apriori  more  about  u,  for  example  that  u  is  analytic  on 
Q  then  we  have 

(1.3)  $(u,?,q,N)  s 


where  C  is  independent  of  N. 

Assume  now  that  we  will  consider  the  family  of  quasiuniform  meshes  as  before 
with  the  uniform  degrees  of  the  elements  but  assume  only  that  pal.  Then  we 
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have  for  u  analytic  on  Q 

- 

(1.4)  $(u.^.q.N)  s  Ce  ^ 

where  C  and  y  are  independent  of  N.  (The  mesh  achieving  (1.4)  is  a 
fixed  one  and  pO')— >0  as  N — >».  ) 

We  see  that  the  rate  of  convergence  depends  on 

a)  properties  of  the  exact  solution  u 

b)  the  admissible  set  of  spaces  S(!IH,q),  OH  6 

We  do  not  know  the  exact  solution  in  advance,  but  for  a  large  class  of 
engineering  problems  we  can  find  a  characterization  of  the  solution  which  can 
be  exploited.  Obviously  the  assumption  that  u  €  H^(£2)  or  u  is  analytic  on 
n  is  an  example  of  a  characterization  of  solutions  of  interest. 

In  practice  we  have  to  consider  a  family  of  problems  characterized  by 
the  useful  data.  Typically  in  the  structural  mechanics  we  deal  with  the 
problem  characterized  by  piecewise  analytic  data,  i.e.  the  boundary  of  the 
domain  is  piecewise  analytic  as  well  the  material  properties,  loads,  etc. 

For  this  class  of  problems  the  assumption  that  u  e  H^(fi)  or  u  is  analytic 
on  £2  is  very  inappropriate.  If  the  domain  £2  has  corners  in  two  dimensions 
or  vertices  and  edges  in  3  dimensions  then  u  is  not  analytic  on  £2  and  u  € 
only  for  small  s  and  the  fact  that  the  solution  u  is  analytic  inside 
Q  is  not  employed. 

In  this  paper  we  will  survey  the  major  results  concerning  the 
performance  of  the  h-p  version  for  solving  elliptic  problems  with  piecewise 
analytic  input  data.  We  will  show  that  for  all  these  problems  the  h-p 
version  leads  to  exponential  convergence  with  respect  to  the  number  of 
degrees  of  freedom.  More  precisely  we  will  show  that  for  a  family  ^  of 
meshes  and  degrees  distribution  q(!in)  we  have 
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(1.5) 


I  u  -  I  „i 

A  AA 


where  a  =  1/3  in  two  dimensions  and  a  =  1/5  in  three  dimensional  problems. 
Constants  C  and  y  >  0  are  independent  of  N.  They  depend  on  the  solution, 
the  domain  £2,  the  distortion  of  the  used  elements,  the  family  f  of  used 
meshes  as  well  the  distribution  qCOD)  of  the  degree  of  elements.  The  family 
of  meshes  which  lead  to  (1.5)  is  a  special  one  and  the  mesh  includes  in  3 
dimensions  the  "needle"  elements  along  the  edges  of  the  domain. 

We  will  discuss  separately  the  two  and  three  dimensional  case  and 
present  the  major  ideas  and  results,  but  without  detailed  proofs  which  will  be 
referred  to.  The  two  dimensional  case  will  be  discussed  in  more  details 
because  it  offers  some  analogies  for  the  3  dimensional  case.  The  numerical 
examples  in  2  and  3  dimensions  will  be  given  also. 

The  survey  papers  about  various  aspects  of  the  h,p  and  h-p  version 
we  refer  to  [1]  [2]  [3]. 
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2.  The  h-p  FE  version  in  two  dimensions 


2.1.  Preliminaries. 


Let  Q  c  =  {x-,x  )  =  {x  |  x.  €  R,  i  =  1,2>  be  a  bounded  domain  with 
12  1 

w  ^  *  N  f  i ) 

the  boundary  Sn.  We  will  assume  that  dQ  =  U  =  F  with  F  being 


Jordan  curves  and 


M(i) 

fd)  .  u  pCi) 

J'l  ^ 


i=l 


i  = 


where  F^.^^  are  analytic  arcs  i.e.  F^.^^  =  {x.  =  <p.  ,(?).  - 

J  J  J-  1  >  J  ^ 

ip.  .(C)  I  c  €  T  =  [-1,  1]>  with  <p.  .(C).  0.  .(?)  analytic  functions  on  T 
l.J  i.J  ^>J 

and  1^  +  1^  'Ai  j(C)|^  ^  a  >  0}.  Here  (^p^  j,  0^  ^)  are 


Cartesian  coordinates  of  T 


By  we  denote  the  open  arc,  i.e.  the 

J  J 


image  of  I  =(-1,1).  The  domain  n  is  R-connected.  We  will  assume  that  the 
1 

curve  F  is  the  boundary  of  the  infinite  component  of  the  complement  of  £2. 

An  example  of  a  domain  Q  of  interest  is  given  in  Fig.  2.1.1. 

Orientation  of  the  curves  is  shown  in  the  Figure  2.1.1  also.  The  endpoints  of 

the  arc  F^.^\  i  =  1,...,R,  J  =  1 . M(i)  are  denoted  by  A^.^^(i.e. 

J  J 

*5-1  =  *5“  =  ‘''i.j''’’  *5” '  4u)’  ““ 

called  vertices  of  £2.  The  measure  of  the  internal  angle  at  Aj^^  is  denoted 

by  We  will  assume  that  0  <  2ir. 

J  J 

We  will  also  understand  the  Jordan  curve  in  an  obviously  generalized 

sense  when  parts  of  different  arcs  could  coincide  so  that  a  slit  domain  could 

•(i) 

also  be  considered  (with  w.  =  Zn) 

J 
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Fig.  2.1.1.  The  scheme  of  a  domain  with  a  piecewise 
analytic  boundary 

Note  that  we  can  have  =  ir  at  some  vertex  (for  example 

1  J  6 

The  vertex  has  to  be  introduced  when  ^j+]  analytic  continuation 

of  Fj^^.  We  can  of  course  place  an  additional  vertex  at  any  place  of 
This  has  to  be  done  because  below  we  will  assume  that  boundary  condition  is 

analytic  and  is  of  the  same  type  on  every  arc. 

TP  4-u  ^  4-  T  II  f^(i)  N_  —  D=  I  I  «(i) 

Further  let  r  =  F  u  F,  F  =  U.  .^F.  ,  F  =  F“  F  =  U.  .^F. 

i,je2)  j  ’  i,J€N  j 

where  D,N  are  the  subsets  of  the  set  {i,j  I  i  =  1,...,R,  j  =  l,...,M(i)}. 

For  simplicity  we  will  assume  that  ^F  ^  0.  Otherwise  we  have  to  add  the 

usual  conditions  of  solvability  and  uniqueness. 

In  Figure  2.1.1  we  have  R  =  2.  We  will  use  R  =  1  in  what  it  follows 

and  will  write  M  instead  of  M(i),  A.  instead  of  A^,^\  etc.  We  note  that 

J  J 
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all  our  arguments  and  results  are  valid  in  the  general  case  too. 

If  the  arcs  are  the  straight  lines  then  we  say  that  £2  is  a  polygon. 

Otherwise  we  speak  about  a  curvilinear  polygon. 

Let  further  be  the  balls  with  the  radius  r .  and  the  center  in 

J  J 

the  vertex  A.  and  let  =  B*!-*  n  £2.  We  will  assume  that  r  <  p,  J  = 

J  J  J  J 

and  p  is  sufficiently  small  so  that  n  F  c  I'j* 

n  =  0  for  every  i  j  and  where  ^  is  the  boundary  of 

J  1  J 

£2J^J^  Finally  we  denote  £2^’’^  =  £2  -  U  with  r  =  (r^, .  .  .  , r^^) . 

The  typical  example  is  shown  in  the  Figure  2.1.2. 


Figure  2.1.2.  The  partition  of  the  domain  £2. 


1  11 
By  H  (£2)  we  denote  the  usual  Sobolev  space  and  Hq(£2)  =  {u  e  H  (£2)  | 

u=0  on  ^F}.  Let  us  now  define  the  space  £0)  c  Hp(£2),  p  =  0  < 

<  1,  i  =  1,...,M.  For  a  =  ^  0  integers,  i  =  1,2,  |a|  =  + 

“2 
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d“u  = 


al^lu 

Xl  X2 


and  for  m  a:  0  an  integer  let 


“u|^(x)  =  ^  |d“u|^(x). 


|a|  =  m 


Further  denote 


$.(x)  =  |x  -  A. I  =  dist  (x,A.),  j  = 

J  J  J 

Let  u  be  such  that  there  exist  constants  C.,  d  ,  and  0  <  p.  <  1, 

J  vJ  J 

j  =  1, . . . ,M,  J  =  0 . M  so  that  for  any  a  =  (a^.a^),  £  0  integers 

i  =  1.2: 


i)  for  X  €  Jlj 


(1.1) 


|d“u1(x)  £  C.d“  a!  l)(x) 

J  J  J 


ii)  for  X  e  Q 


(r/2) 


(1.2) 


|d“u| (x)  s  Cgd®  a! 


where  d.  =  (dl^^  d^.^^  >  1.  d“  =  (d^^  ^  )“^  (d^^^  )“2,  a!  .=  a  !a  !.  0!  =  1. 

J  J  J  v3  s3  J 


By  SO)  we  denote  the  set  of  all  functions  u  satisfying  (1.1),  (1.2). 

Let  us  underline  that  overlaps  the  domains  Hence  in 

this  overlap  area  the  function  u  €  SO)  under  consideration  satisfies  both 

conditions  i)  and  ii).  We  always  will  assume  an  overlap  without  mentioning 

( r ) 

it  explicitly  and  write  often  only.  The  description  of  the  space 

SO)  depends  on  r,  nevertheless  this  is  not  essential  because  the  constants 
C  and  d  in  (1.1)  and  (1.2)  are  not  explicitly  specified.  On  the  other  hand 


the  dependence  on  /3  is  essential.  The  assumption  0  <  <  1  guarantees 


that  SO)  c  H^(n)  and  SO)  c  C°(n). 
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2.2.  The  boundary  value  problem. 

Let  us  consider  the  problem 
(2. la)  -A  u  =  f  in  n, 

(2. lb)  u  =  0  on  ^r, 


(2.1c) 


3u  N„ 

=  g  on  r. 


We  understand  the  problem  (2.1)  in  the  weak  sense.  Find  u  e  such  that 

1 

for  any  v  €  H^CQ) 


B(u,v) 


r  du  dv  ^  du  dv 


holds. 

We  will  assume  that  the  functions  f  and  g  are  such  that  the  solution 
u  €  £(p)  for  some  p  dependent  on  Q,  more  precisely  depends  only  on 

the  internal  angle 

In  [4]  [5]  [6]  we  have  analyzed  in  detail  the  space  of  input  data  g 
and  f  for  which  u  €  £(p)  and  gave  a  detailed  description  of  these  spaces. 
As  special  case  we  have  shown  that  u  €  £0)  with  properly  selected 
dependent  on  if 

a)  function  f  is  analytic  on  Q 

-  N- 

b)  g  is  analytic  on  every  arc  F.  c  F 

J 


Remark  2.1.  In  [4]  we  have  analyzed  the  regularity  of  the  solution  u  of  the 
boundary  value  problem  for  general  elliptic  operator  with  analytic  co¬ 
efficients  on  Q  and  have  described  the  spaces  of  g  and  f  that  u  €  £0). 
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Remark  2.2.  In  [6]  we  have  analyzed  the  regularity  of  the  solution  u  of  the 
elasticity  problem  and  have  shown  that  u  €  SO)  when  input  data  are 
piecewise  analytic. 

Remark  2.3.  In  [7]  we  analyzed  the  eigenvalue  problem.  We  have  shown  that 
the  eigenfunctions  belong  to  the  space  £(p). 

Remark  2.4.  We  restricted  ourselves  to  the  case  that  u  =  0  on  In  the 

above  mentioned  papers  we  have  shown  (as  special  case)  that  if  u  =  g 
with  g  analytic  on  the  arcs  Fj  €  ^F  we  get  also  u  €  SO).  The 

restriction  u  =  0  on  ^F  will  simplify  our  approximation  theorem.  In 

general  when  u  ^  0  on  ^F  we  have  to  replace  the  nonhomogeneous  Dirichlet 
boundary  condition  by  a  properly  selected  sequence  of  functions.  We  will  not 
address  here  this  problem.  (See  for  more  e.g.  in  [8]). 

Let  us  underline  that  in  the  engineering  practice  we  essentially  always 
deal  with  piecewise  analytical  data  and  hence  the  solution  under  consideration 
always  belongs  to  the  space  SO). 

It  is  essential  to  deal  with  the  spaces  of  solutions  which  are  as  small 

as  possible  but  will  include  practically  all  the  cases  of  engineering 

importance. 

2.3.  The  elements 

As  usually,  we  will  assume  that  the  domain  Q  is  covered  by  the  mesh  JR 
which  partitions  the  domain  Q  into  elements  J  and  will  consider  a  family 
^  of  admissible  meshes  JR  e  3^.  Hence  we  will  write  3R  =  {S'}.  We  will 
consider  first  family  9  of  meshes  with  quadrilateral  elements. 

2 

Let  us  denote  by  D  the  standard  quadrilateral  element,  D  =  I  ,  I  = 
(-1,1)  i.e. 
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D  =  t)2  1  hjl  <  1>  h2l 

We  will  assume  that  any  J  €  OH  is  an  image  of  D  by  the  one  to  one  mapping 


My,  i.e. 

y  =  I  T}^,  T}^  e  D>. 

About  i  =  1,2  we  will  assume  that  there  is  viV^)  >  0  such  that 

i)  for  any  |a|  =  m  >  0  and  any  t)  e  D,  f  =  1,2 

(3.1)  |d“x^(t))|  £  Cgd™  m!  i^,  m  =  1,2... 


ii)  The  Jacobian  determinant  7  satisfies 


(3.2) 


^  1^1 


d(X^,X^) 

dlTj^Tr?^ 


2 


In  (3.1)  and  (3.2)  v  dependents  on  the  element  J  but  the  constants 
CQ,Cj,Cj^,dQ  are  indepndent  of  S’  €  3)1  €  Obviously  u  is  related  to  the 
size  of  the  element  S’. 


Remark  3.1.  If  the  boundary  of  the  element  lies  on  the  boundary  F  then 
there  is  obviously  a  relation  between  the  constants  in  (3.1)  (3.2)  and 
description  of  the  analytic  arc  F.. 

Remark  3.2.  We  will  assume  that  every  vertex  of  Q  is  also  the  vertex  of  an 
element. 

Remark  3.3.  Let  us  note  that  from  (3.1)  we  see  that  XAi\)yi  =  1,2  are 
analytic  functions  on  D  and  hence  they  can  be  extended  beyond  D  in  some  of 
neighborhood  of  D, 

Remark  3.4.  Our  assumptions  imposed  on  the  elements  guarantee  that  the 
aspect  ratio  of  S'  €  OH  e  S'  is  bounded. 
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Obviously  any  quadrilateral  with  the  straight  sides  satisfies  the 
conditions  mentioned  above.  For  example  for  i  =  1,2  of  order 

one  and 

Xi  =  V’(a^(l+r)j)  (1-7/2) 

+  bj^  (1+7/^)  (1+7/2)  +  c^(l-7/j)  (1+7/2 ) ) 

maps  D  onto  7  (shown  in  the  Figure  2.3.1)  provided  some  conditions  on  the 
coefficients  a^^.b^.c^  are  imposed. 


Fig.  2.3.1.  The  element  S' 

Let  now  S'  be  a  triangular  element.  Then  we  will  assume  that  7  is 
the  image  of  the  tnaster  element  T  =  •(t/^  ,7/2 1  (t/^  ,7/2)  e  D,  ■'>2  ”'^1^  c  D  by  the 

mapping  which  satisfies  (3.1)  (3.2)  and  ^^.(D)  c  0.  Hence  S  is  a 

half  of  My(D). 

We  now  describe  the  elements  of  the  meshes  we  will  consider. 

1)  Elements  in 

j 

We  will  consider  two  kinds  of  elements  7  and  S'q  separately. 

The  element  .  Let  J  c  then  A.  S.  Denote 

J  J 

(3.3a)  ic.(y)  =  min  $.(x), 

xeJ  ^ 
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(3.3b) 


with  C,  C  independent  of  9"  e  3H  e 

Condition  (3.8)  shows  that  there  is  only  finite  number  of  elements  in 
(dependent  on  various  constants  in  (3.1)  (3.2),  (3.8)  but  independent 
of  jn  €  S').  Elements  of  3)1  are  curvilinear  with  the  sizes  proportional  to 
the  distance  from  the  vertex  (except  those  elements  which  contain  the  vertex). 


2.4.  The  meshes  and  the  finite  elements. 

Let  us  consider  a  family  ^  of  meshes  on  £2.  A  mesh  !II1  e  3"  is  a 

partition  of  £2  into  the  set  of  elements  3"  €  OH  which  satisfy  conditions 

given  in  the  Section  2.3  and  £2  =  U  To  every  S’  €  OH  we  have  associated 

JeOH 

the  analytic  map  M^.  Obviously  we  can  speak  about  the  vertices  and  sides  of 

3".  We  will  assume  that  the  elements  of  3"  €  OH  satisfy  the  usual  conditions 

i.e.  if  S’j,  3’j  €  OH  then  either  are  disjoint  or  have  common  vertex 

or  common  sides.  If  W.,  W.  have  common  sides  then  the  mapping  ,  Mor 

1  J  -J  i  J  i 

coincides  on  the  common  sicje. 

( r  * ) 

Let  us  now  consider  the  mesh  of  elements  In  ^  ,  Because  (3.3) 
and  (3.4)  we  can  speak  about  the  layers  of  the  elements.  The  zero  layer 


consists  of  all  elements  which  (closure)  includes  the  vertex  Aj.  Denote  by 
the  k-th  layer  of  elements.  Then  we  define  layer  as  the  set  of  all 


elements  S' 


c  £2^.-J^ 


such  that  they  do  not  belong  to  the  layers 


J  = 


1,2 . k-1  but  3"  A  U  3’  ^0.  From  the  conditions  (3.1)-{3.5)  we  see 

that  there  is  a  constant  dependent  only  on  the  constants  in  (3.1)-(3.5) 

so  that  any  layer  !£.^  consists  of  at  most  Kq  elements. 

We  will  now  consider  the  family  3-  of  meshes  (of  the  elements  described 
in  Section  3.2)  which  is  characterized  by  the  size  of  the  smallest  element 
in  the  mesh.  We  will  assume  for  given  <r  <  1  and  n  >  0  an  integer  we  have 
for  c  £2^.’'j\  A.  e 


(4.1) 


C  O'"  s  k(3'q'^^  s  C 


Hence  the  mesh  3H  is  characterized  by  the  parameters  (cr.n)  and  we 


will  write  !)H((r,n).  Note  that  with  increasing  n,  the  number  of  elements  in 

£2^.’'J^  is  growing,  while  the  number  of  elements  in  is  uniformly 

J 

bounded.  In  fact  we  may  assume  that  the  mesh  in  is  independent  of 
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n.  We  show  the  typical  mesh  in  the  neighborhood  of  a  vertex  in  the 

Fig,  2.4.1.  The  layers  of -the  elements  are  shadowed.  In  general  ((r,n)  can 

depend  on  j,  i.e,  can  be  different  in  every  vertex  neighborhood. 

Nevertheless  we  will  assume  that  ((r,n)  are  independent  of  .  j. 


Fig.  2.4.1.  Typical  mesh  with  shadowed  layers 

Mesh  of  this  type  will  be  called  a  geometrical  mesh  with  the  factor  <r.  We 
can  also  have  cr  =  1.  Then  the  mesh  JIKor.n)  will  have  number  of  elements 
independent  of  n.  We  can  then  assume  that  the  mesh  is  a  fixed  one.  If  we 
know  apriori  that  in  the  neighborhood  of  some  vertices  no  singularity  will 
occur,  then  we  select  cr  =  1  in  these  neighborhoods.  This  is  typical  for 
example  when  we  impose  the  symmetry  condition  on  some  part  of  the  boundary. 

Remark  4.1.  Denote  by  the  element  in  the  k  layer.  Then  using  (3.1) 

(3.3) (3.4) (3.5)  and  (4.1)  we  have 

(4.2)  <  C<r^p"^ 

with  <  1,  t  -  1,2  (depending  on  the  constant  in  (3.1)-(3.5).  In  practice 
we  mostly  construct  the  meshes  such  that  (3.1)-(3.5)  hold  together  with 
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(4.3) 


i . e .  we  have 
(4.4) 


Pj  =  O',  i  =  1,2 


_ ( k )  ( k ) 

Remark  4.2.  Denote  by  k  resp  k  the  maximum  resp  minimum  of  icO") 

over  all  elements  in  the  layer.  Then  s  r.  and  hence  number  of 

^-k  J 

layers  does  not  exceed  Cn  whether  (4.1)  (resp  (4.2))  or  (4.3)  is  used. 

The  shape  functions  of  degree  p(J)  are  as  usually  given  on  the  master 
element  D  or  T.  On  D  we  will  assume  that  the  degree  is  separately  in  every 
variable.  For  simplicity  we  will  assume  that  the  degrees  of  the  elements  are 
uniform  although  the  nonuniform  distribution  is  more  effective,  but  would  have 
the  Soime  rate  of  convergences. 

Finally  we  will  denote  by  3(311, p)  =  S(3R(CT‘,n),  p)  =  S(<r,n,p)  c  Hq(Q)  the 
finite  element  space  under  consideration. 

2.5.  The  h-p  version  of  the  finite  element  method 

The  finite  element  method  for  the  model  problem  (2.1)  reads: 


Find  Ug  e  S  =  S(3R,p)  =  S(<r,n,p)  such  that 


(5.1) 

where 

(5.2b) 

(5.2b) 


B(Ug,v)  =  F(v),  Vv  e  S((r,n,p), 


B(Ug,v)  = 


F(u)  = 


fdus  dv  ^  8us  dv 
[dxi  dxi  5x2  9x2^ 


f  V  + 


g  V. 


N, 


Then  with 


(B(u,u))^''^  =  Hull 
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we  have 


>1^0  ■  ^s“e  =  "^0  ■  ^"e 

;«eS 

where  is  the  exact  solution  defined  in  the  Section  2.2.  Obviously 

1 

||u|L  =  luL  where  by  |u|,  we  denoted  the  seminorm  in  H  (0). 

Hi  X  X 

Let  us  now  formulate  typical  theorems.  We  outline  the  main  ideas  of 
their  proof  in  the  Section  2.6.  For  details  we  refer  to  [9],  [10],  [11]. 


Theorem  2.5.1.  Let  u^  be  solution  of  the  problem  (2.1).  Assume  further 
that  for  any  x  e  £2  and  any  a  =  (a^.a^),  i  0  integers  we  have 
(see  (1.2)) 

(5.3) 

Consider  now  the  mesh  !IIl(<r,n)  with  cr  =  1  and  p(y)  =  n,  9"  e  3H.  Then 


o 

Q 

<  C  d^ 

with 

(T  =  1 

a 


(5.4a) 


(5.4b) 


A;€S(o-,n,p) 


s  C  e 


-rn 


ajeS((r,n,p) 


-rN 


1/2 


where  T  >  0  and  C  >  0  depend  on  the  solution,  the  elements,  the  domain 
but  are  independent  of  n  and  N;  N  is  the  number  of  degrees  of  freedom 
(N  =  dim  S(a',n,p) ) . 


Remark  5.1.  The  condition  (5.3)  is  equivalent  with  the  assumption  that  u  is 
analytic  on  Q, 


Remark  5.2.  Because  (t  =  1  the  number  of  elements  is  independent  of  n. 
Hence  the  method  is  the  p-version  of  the  finite  element  method. 


16 


Theorem  2.5.2.  Let  be  the  solution  of  the  problem  (2.1).  Assume  that 

Uq  €  SO)  as  defined  in  the  Section  2.1.  Consider  now  the  geometric  mesh 
!)Il(<r,n)  with  a-  <  1  and  pCJ)  =  pn,  p  >  0,  or  properly  chosen.  Then  there 
exist  C,  y,  3r  >  0  such  that 


(5.5a) 


;t;€S((J',n,p) 


(5.5b) 


l^o-^slHi(n) 


inf  |Uq 

;»;€S((r,n,p) 


£  Ce 


where  C,y,  y  are  independent  of  n  and  N. 


Remark  5.3.  If  (t  <  1  then  number  of  elements  in  the  mesh  3H(<r,n)  is 
proportional  to  n,  and  for  p  =  pn  (5.5b)  follows  directly  from  (5.5a). 

In  the  Theorem  5.2  the  assumption  u^  €  SO)  is  essential.  As  we 
mentioned  in  the  Section  1,  this  assumption  is  satisfied  practically  for  every 
problem  in  structural  mechanics  where  all  the  input  data  are  piecewise 
analytic. 

The  h-p  version  simultaneously  changes  the  mesh  and  the  degree  of 
elements.  In  practice  often  the  geometrical  mesh  JIKo'.n)  is  a  priori 
constructed  for  fixed  o*  and  n  =  n^  and  then  p  is  increased  i.e.  the 
space  S((r,nQ.,p)  with  p  increasing  is  used.  Then  we  see  two  phases  in  the 
behavior  of  the  method.  In  the  first  phase  the  convergence  is  very  similar  as 
in  the  h-p  version  with  the  exponential  convergence  and  in  the  second  phase 
for  large  p  the  performance  is  similar  as  in  the  p-version  with  an 
algebraic  rate  of  convergence. 

Remark  5.4.  The  constants  y  depend  among  others  on  <r  and  p.  In 
[12]  we  analyzed  one  dimensional  case  and  found  that  <r  »  0.15  is  optimal. 

In  the  next  section  we  will  see  some  numerical  results  related  to  the  question 


of  the  optimal  selection  of  c. 

Remark  5.5.  The  h-p  version  for  elliptic  equations  of  order  2m  was 
addressed  in  [13]. 

Remark  5.6.  We  formulated  theorems  1  and  2  for  Laplace  equation  only. 
Because  the  essential  feature  of  these  theorems  is  the  approximation,  they 
hold  in  general  setting  when  u^  €  £0),  for  example  for  the  elasticity 
problem. 

2.6.  Numerical  experiments 

In  this  section  we  will  show  a  typical  numerical  example.  Consider  the 
elasticity  problem  on  a  cracked  domain  shown  in  the  Fig.  2.6.1a.  Because  we 
will  consider  the  symmetric  problem,  only  the  half  domain  will  be  considered 
as  show  in  Fig.  2.6.1b. 


Fig.  2.6.1.  The  scheme  of  the  cracked  domain 
a)  The  cracked  domain,  b)  the  symmetric  half  of  the  cracked  domain. 
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We  will  assume  that  material  is  isotropic  homogeneous  with  Poisson  ratio  v  = 

0.3.  On  the  boundary  F  of  Q  depicted  in  the  Fig.  2,6.1,  we  impose 

nonhomogeneous  traction  (Neumann)  conditions  so  that  the  exact  solution  is 

“1/2 

symmetric  stress  intensity  mode  for  which  the  stress  behaves  as  0(r  ). 

Hence  we  have  only  singularity  in  the  origin  in  the  vertex  A^.  Therefore  the 
mesh  will  be  refined  only  in  the  neighborhood  of  the  vertex  (i.e.  <r  =  1 
in  the  neighborhood  of  all  other  vertices).  Because  of  the  symmetry  we  will 
consider  only  the  problem  on  the  half  of  the  domain  as  shown  in  Fig.  2.6.1b. 

We  will  use  o'  =  0.15. 

In  the  Fig.  2.6.2  we  show  the  sequence  of  the  meshes  with  n  layers 
(not  drawn  in  scale). 


MESH  Mg  MESH  M3 


Fig.  2.6.2.  The  sequence  of  the  geometric  meshes  JIKcr.n),  n  =  0 . 5. 
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We  will  use  p  =  n  +  1  in  our  example  and  xmiform  p.  Based  on  the 

(n) 

Theorem  2.5.2  we  expect  that  the  error  He  is  bounded  by  C  e 

with  C  and  y  independent  of  n  and  N.  We  will  assume  that 


(6.1) 


|e  II 


ER 


II  (n^i 


C(n)e 


-yN 


1/3 


and  determine  the  value  y  from  the  two  successive  values  of  ||e^^^||^. 

Then  we  compute  for  every  n  the  value  C  and  y.  The  results  are  given  in 
the  Table  2.6.1.  Although  Theorem  2.5.2  gives  only  an  upper  estimate  it  seems 
that  y(n)  and  C(n)  converge  as  n— ^oo.  This  is  not  surprising  in  our 


particular  case,  but  it  cannot  be  concluded  neither  from  the  theorem  or  the 
theory  presented  here. 


Table  2.6.1.  Performance  of  the  h-p  version 


n 

P 

N 

N1/3 

"'""•er’'- 

y(n) 

C(n) 

0 

1 

9 

2.08 

60.92 

0.741 

1.455 

1 

2 

48 

3.63 

20.23 

0.740 

1.455 

2 

3 

121 

4.95 

7.61 

0.776 

2.098 

3 

4 

256 

6.35 

2.57 

0.720 

1.810 

4 

5 

477 

7.82 

0.90 

0.670 

1.683 

5 

6 

808 

9.31 

0.33 

0.670 

1.688 

In 

the  Fig.  2. 

.6.3  we 

(n ) 

show  the  error  |le^  \\rT>  sis 

EjK 

function 

of  N.  We 

plot  the  graph  in 

N  X 

Ig  lle^'^^ii^  ; 

4  scale.  Then 

the  graph  would  be  the 

straight  line  if  the  error  would  obey  exactly  (6.1). 
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N1/3_^ 

Fig.  2.6.3.  Performance  of  the  p  and  h-p  version 
for  S(<r,n,p),  <r  =  0.15 

In  the  Fig.  2.6.3  we  also  show  the  error  of  the  p-version  i.e.  the  error 
using  S(0.15,  n,  p)  for  fixed  n  and  1  <  p  <  8.  The  h-p  version  graph 
connects  the  points  giving  the  accuracy  for  p  =  n  +  1.  We  see  that  we  obtain 
a  straight  line.  Table  2.*6.1  and  Fig.  2.6.3  shows  that  the  (6.1)  gives 
not  only  the  upper  estimates  but  also  describes  well  the  behavior  of  the 
error.  From  the  Fig.  2.6.3  we  also  see  that  for  given  n  the  error  behavior 
as  function  of  p  can  be  divided  into  two  phases.  In  the  first  phase  when  p 
n  +  1  we  see  exponential  convergence  while  for  p  ^  n  +  1  we  see  only  an 
algebraic  one. 

The  performance  depends  among  others  on  the  value  of  cr  and  p.  In  one 
dimensional  setting  o'  «  0.15  is  optimal,  see  [12].  In  the  Fig.  2.6.4  we 
show  the  performance  of  the  h-p  version  for  p  =  n  +  1  as  function  of  cr. 
The  Fig.  2.6.4  is  drawn  in  the  same  scale  as  Fig.  2.6.3.  We  see  that  the 
value  O'  -  0.15  gives  the  optimal  results. 


o 


N 


50 

30 

20 

t 

^  5 

_5 

Q)  2 

1 

0.5 

0.3 


N1/3 — ► 


Fig.  2.6,4.  The  performance  of  the  h-p  version  for  various  cr 


Remark  6.1,  We  have  shown  the  performance  of  the  h-p  version  as  function 


of  N.  Of  course  for  higher  p  the  stiffness  matrix  is  more  dense  and  the 


cost  of  construction  of  the  stiffness  matrix  is  higher  too.  Hence  the  right 


performance  description  will  be  graph  of  the  computer  cost  x  accuracy.  This 


problem  was  addressed  in  [14]  [15]. 


Remark  6.2.  The  meshes  for  the  h-p  version  have  a  special  character.  The 


usual  available  mesh  generators  are  not  producing  the  appropriate  meshes  for 


the  h-p  version  yet. 
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2.7.  Outline  of  the  theory. 

In  this  section  we  will  briefly  outline  the  theory  leading  to  the  theorem 
2.5.1  and  theorem  2.5.2.  The  main  ideas  in  two  dimensional  and  three 
dimensional  case  are  similar,  but  technicalities  are  much  more  difficult  in 
3  dimensions. 

1)  The  regularity  problem. 

The  first  major  problem  is  to  characterize  the  space  of  solutions  under 
consideration.  The  space  should  be,  on  one  hand,  as  small  as  possible,  to 
give  the  possibility  to  employ  its  special  properties,  but  on  the  other 
hand  the  space  has  to  be  large  enough  that  it  would  cover  most  problems  in 
engineering  practice.  As  was  mentioned  earlier  typical  engineering  problems 
are  characterized  by  piecewise-analytic  input  data.  Hence  we  need  to  describe 
the  spaces  of  the  solutions  of  such  problems. 

There  is  a  large  literature  about  the  regularity  of  the  solution  of 
elliptic  problems.  Nevertheless  those  theories  are  not  directed  to  the  goals 
of  numerical  solutions.  In  a  series  of  papers,  [4]  [5]  [6]  [7]  we  developed 
the  theory  which  characterizes  the  regularity  of  the  solutions  in  the  terms  of 
countably  normed  spaces  (see  (1.1),  (1.2))  which  is  very  advantageous  for  the 
analysis  of  the  h~p  version.  This  theory  in  the  cited  papers  encompasses 
larger  class  than  the  one  used  here. 

We  do  partition  the  domain  in  the  (overlapping)  areas  and  characterize 
the  regularity  in  these  areas  in  a  special  way.  In  two  dimensional  settings 
we  distinguish  between  vertex  neighborhood  (denoted  above  by  £2^T^^)  and  the 
internal  domain  (denoted  by  ) .  In  3  dimensions  we  have  considered  four 

regions  separately. 

2 

2)  Approximation  of  functions  defined  on  the  standard  square  D  =  I  ,  I  = 

(-1,1). 
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The  major  theorem  is 

Theorem  2.7.1.  Let  d“u  e  L^(D),  a  =  ,  0  s  <  t^  +  1  ^  0, 

I  =  1,2.  Then  there  exists  a  polynomial  PCx^.x^)  of  degree 
Xj^.x^  such  that  for  0  £  s  1  and  1  ^  s^  s  t^,  £  =  1,2  we  have 


=4TtTTir?zn^  I 

OsSpSl 


(7.1) 


(ta-sa) ! _ 

(ta+sa  +  2(l~<xa))! 


O^pi^l 


For  0  s  ^  1,  =  0 


(7.2)  IID“(u-P)(tl,X2)l|2,,j,  C  (t34^^1(l-«,))T  Z 

0£6i:£l 


llnPi  •  ®2'*’1,,(|2 

IID  uIIl2(D) 


For  0  s  s  1,  =  0 


(7.3)  ||D“(U-P)(x^.±l)llg3„)  »  C  (t,4i'r2(l-«.))T  Z  "‘’°‘*'’^''‘"l=(D) 

L  OsBa^l  J 


(7.4) 


(u  -  P)  (±1.  ±1)  =  0 


tti+aa 

We  denoted  =  — - —  in  (7.1)  -  (7.3)  and  thereafter. 

ax“>  ax«^ 

We  see  that  smoothness  of  the  trace  of  the  function  (u  -  P)  and  its 
norm  on  the  sides  of  D  is  estimated  by  various  derivatives  of  u  in  D.  As 
a  simple  corollary  from  the  theorem  7.1  we  get 

Theorem  2.7.2.  Let  i/»(x)  =  (u  -  P(x^,X2)),  x^  =  X2  =  x  e  I  then  for 
la|  =  1,  (i.e.  X  lies  on  the  diagonal  of  D), 
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HD  -  C  (tj+sj)! 


■(^  I 

0sa2^1 


l2(D) 


(7.5) 


( t2-S2 ) t 
( t2+S2 ) ! 


OSttiSl 


The  constant  C  in  Theorem  2.7.1  and  2.7.2  is  independent  of  t,s.  For  the 
proof  in  the  3  dimensional  setting  we  refer  to  [16]. 

3)  Regularity  of  the  transformed  function  on  the  standard  square 

As  said  above,  there  is  a  mapping  which  maps  D  onto  S'.  We 

assumed  that  the  mapping  has  properties  given  in  (3.1)  (3.2).  We  also 

have  assumed  that  the  solution  u  €  SO).  Consider  now  the  geometrical  mesh 

!)R(<r,n)  and  let  7  €  and  W  r\  k  =  0.  Then  for  x  e  7  we  have  from 

J  J 


(1.1) 


|D%|(X)  = 


(7.6) 


=  a!  k'-^, 

J  J 


where  we  wrote  k  instead  of  Kj  defined  in  (3.3),  dj,d2  instead 

d^.^\  d^.^\  and  p,  C  instead  of  p.,  C.. 

J  J  J  J 


Let  now 


(7.7) 


LI(7))  =  u(Mg.(T))  )  , 


then  U(i))  is  defined  on  D  and  we  have 


Theorem  2.7.3.  For  any  integer  s  s  0  and  any  tj  €  D  we  have 


For  0  £  s  1 


(7.8) 


|D®’“2u(t7)1  £  Cd®s! 


For  0  £  £  1 
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(7.9) 


U(d)  I  s  Cd®  s! 


where  the  constants  C,d  are  independent  of  v  and  s. 

Quite  similar  result  holds  for  J  The  elements  for  which 

n  Aj  0  have  to  be  treated  separately  as  shown  below. 

4)  The  approximation  on  the  element  3^ 

We  use  now  the  estimates  (7.8),  (7.9)  together  with  the  results  of  the 
Theorem  2.7.1  and  obtain 

Theorem  2.7.4.  There  exists  a  polynomial  P(Vj^,'n2^  degree  p  in 

and  such  that  for  Osa^sl,  i  =  l,2 


(7.9) 


||d“(U-P)||^2(d)  F(d)V"^^^  . 

£=1 


(7.10a) 


HD' 


1^2 


(U-P)  (±1.7,2)11^2(1)  =£  ^^F(d))V  , 


(7.10b) 


I1d“^  (U 


-  P)(7,j.  ±1)11^2(1)  ^  ^^(F(d))V  . 


(7.11)  (U  -  P)(±l.  ±1)  =  0  . 

Function  F(d)  is  increasing  but  for  any  d,  F(d)  <  1.  The  function  F(d) 
is  obtained  from  (7.1)  by  the  optimal  choice  of  s  for  given  degree  of  the 
polynomial.  For  detail  properties  of  this  function  we  refer  to  [9],  [10]. 
From  (7.9)-(7.11)  we  see  that  the  convergence  rate  in  p  is  exponential. 
Using  now  the  transformation  and  utilizing  (3.2)  we  get 

-1 

Theorem  2.7.5.  Let  y  =  P,  where  P  is  the  polynomial  in  the 
Theorem  2.7.4.  Then 
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i  =  1.2 


^  Ia|=l 

(7.12) 

<  Ck  2(1-P)qP 

where  0<Q<1.  It  is  independent  on  k,  but  depends  on  the  constants  in 
(3.1)  and  (3.2) 

The  elements  which  have  a  vertex  in  the  A .  have  to  be  treated 

J 

separately.  Using  the  assumed  properties  of  u  it  is  easy  to  obtain 


Theorem  2.7.6.  Let  €  JR  is  such  that  n  Aj  ^  0. 

(7.13)  "as;  ‘'^-'‘>“^(70)  = *  = 

Function  <p  is  here  mapped  bilinear  function. 


Then  we  have 

1,2 


5)  The  adjustments  of  local  approximation 

The  results  in  4)  allqw  to  estimate  the  approximation  by  a  polynomial 
in  every  element  separately.  Nevertheless  this  construction  of  the 
approximation  does  not  lead  to  the  function  which  is  continuous  across  the 
boundaries  of  the  elements  i.e.  conforming  elements.  The  continuity  is 
guaranteed  only  in  the  vertices  of  the  elements  because  the  approximation 
coincides  here  with  the  approximated  function.  Hence  we  have  to  make 
corrections  which  will  delete  these  (sides)  discontinuities  .  This  is  made 
as  follows.  Let  and  elements  with  common  edge  I  and  let 

£  r>  r  =  0.  Assume  that  and  are  the  sides  of  the  standard  square  D 

which  by  Ms-  and  Mo-  are  mapped  on  the  edge  1. 

Ji  J2 

Without  any  loss  of  generality  we  will  assume  that  =  S.  Let 

Ui*  i  =  1,2  be  the  mapped  functions  U  and  the  polynomial  P  considered 

in  the  Theorem  2.7.4.  Then  obviously  “  ^2  ^  0  =  Pj  -  P2  is  a 
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polynomial  of  degree  p  on  S  and  0=0  in  the  end  points  of  S.  Using 
(7.10  a)  we  get  with  0  <  Q  <  1 


(7.14) 


H  (S)  * 


Without  any  loss  of  generality  we  will  assume  a  and  construct  on  D  a 
polynomial  of  degree  p  (in  t)j  and  tj^)  such  that  =  0  on  S  and 
Qj  =  0  on  SD  -  S  and 


(7.15) 


HQ.  II  1  s  CII0II  ,  . 

H  (D)  H  (S) 


We  use  now  instead  the  polynomial  Pj“Qj  8©^ 

(7.16)  IIU  -(P,-Q,)||  ,  £  IIU  -P  II  +  HQ  II 

^  ^  ^  r(D)  ^  ^  H^(D)  ^  H^(D) 

i  C(HU  -PJI  .  +  II0H  .  ) 

H  (D)  H  (S) 

.  C(.f a")'''". 


Hence  the  error  of  the  finite  element  method  can  be  estimated  by  the  errors  in 
the  single  elements  as  described  in  Theorem  7.5  and  7.6. 


6)  Proof  of  the  Theorem  2.5.1,  and  2.5.2. 

Let  us  first  address  Theorem  2.5.1.  By  the  assumption  the  mesh  III  has 
m((r,n)  elements  (independent  of  m)  and  u  satisfies  (5.3).  Using  the 
results  in  the  previous  paragraphs  we  construct  function  (p  e  S((r,n,p),  cr  =  1 
such  that 

llu-»)||^  ^  CmZ*^  =  Cme”^“P,  a  >  0,  0  <  Z  <  1, 

which  yields  (5.4a).  Because  the  number  of  degrees  of  freedom  N  (=  dimension 

2 

of  S(<r,n,p))  is  of  order  mp  we  obtain  (5.4b). 


Let  us  now  address  Theorem  2.5.2.  We  divide  the  elements  of  the 


geometrical  mesh  in  three  groups.  Elements  in  elements  in 
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J  =  1 . M  but  such  that  J  n  A.  =0  and  J  c  S'  r\  A  *  a.  Obviously 

J  J  J 

we  have  at  most  CM  elements  (with  M  being  the  number  of  vertices)  of  the 

third  group.  Using  Theorem  2.7.6  the  error  (square)  on  these  elements  in  the 

third  group  is  estimated  by  CM((r^)^^^  where  p  =  max  In  the  second 

group  in  every  we  have  at  most  Cn  layers  and  the  number  of  elements 

J 

in  every  layer  is  uniformly  bounded  (see  Sec.  2).  Using  Theorem  2.7.5  we  see 
that  the  error  (square)  in  the  domain  ^  bounded  by 

*^^^2n(l-3j)^-2(l-p)jQp  ^ 

j=l 

^  ^,^2n(l-p)  p-2n(l-p) 

where  Q  <  1  is  independent  of  p  and  o'.  Hence  the  total  error  (square)  in 
the  entire  Q  (composed  from  the  errors  of  the  three  mentioned  graphs)  is 
for  p  =  p  n  bounded  by 

^  ^^)  ]  ■  ^  ^ 

provided  that  p  was  properly  chosen  depending  on  <r,  p,  p. 

Obviously  the  number  of  elements  in  3H  is  bounded  by  Cn  and  we  get 

(5.5a),  (5.5b)  then  follows  from  the  fact  that  the  dimension  of  the 

„  .  2 

space  of  polynomials  of  degree  pm  u  is  p  . 

So  far  we  have  assumed  that  the  elements  are  quadrilateral.  If  they  are 

triangular  then  we  use  first  the  approximation  on  the  square  D  and  the  only 

difference  is  in  the  adjusting  phase.  Here  we  use  Theorem  2.7.2  for  the 
estimate  on  the  diagonal  of  the  square  D.  The  extension  from  the  sides  of 
the  master  triangle  inside  it  is  standard.  Of  course  we  have  to  see  that  at 

the  diagonal  of  D  the  function  is  the  polynomial  of  degree  2p,  but  this 

influences  the  constants  but  not  the  rate  in  the  bounds. 
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In  this  section  we  will  formulate  the  results  for  the  3  dimensional 
problem  which  are  analogous  to  those  mentioned  in  Section  2.  We  will 
underline  similarities  and  differences  between  the  two  and  three  dimensional 
results. 

3. 1.  Preliminaries 

Because  of  simplicity  we  will  consider  only  polyhedral  domains  although 
our  results  have  much  more  general  character. 

Figure  3.1.1  shows  a  typical  domain.  The  planes  KLNO  and  MLNO  coincide 
(crack  type).  By  A,B...  we  denoted  the  vertices  of  the  domain.  The  edges 
of  the  domain  Q  are  straight  lines  with  the  vertices  at  their  end.  For 
example  AB  is  an  edge.  We  have  in  the  Figure  3.1.1  only  edges  where  the 
integral  angle  ^  tt.  Nevertheless  for  the  same  reasons  as  in  2  dimensional  case 
we  can  or  have  to  add  additional  vertices  or  edges,  e.g.  the  edge  BL  etc.  A 
set  of  edges  creates  the  boundary  of  a  face,  which  is  a  polygon.  For  example 
BCHI  JKLMG  are  the  vertices  of  a  face.  As  in  the  two  dimensional  case 
we  will  assume  that  on  every  face  a  boundary  condition  is  of  the  same  type. 

In  our  case  we  have  only  plane  faces.  Nevertheless  the  edges  and  faces  can  be 
curved  too. 


P  Q 

Fig.  3.1.1.  The  polyhedral  domain 
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In  the  two  dimensional  case — in  the  Section  2.1 — we  partitioned  the  domain  in 
the  neighborhood  of  the  vertices  and  rest  of  the  domain.  The  three 
dimensional  case  is  more  complex.  Here  we  have  to  partition  the  domain  in  the 
following  types  of  regions. 

1)  Neighborhood  of  the  edges  (not  close  to  the  vertices). 

2)  Neighborhood  of  the  vertices  edges  (close  to  the  vertices  and  edges). 

3)  Neighborhood  of  the  vertices  (not  close  to  the  edges). 

4)  The  regular  region. 

Let  us  describe  these  regions  more  precisely. 

R  5 

1)  The  neighborhood  of  the  edge  Q^’  (not  close  to  a  vertex). 

Assume  that  the  edge  e  under  consideration  is 

(1.1)  «  =  {x^,X2,x^  I  x^  =  x^  =  0,  0  <  x^  <  1>. 

Then  we  denote  forR>0,6>0 

(1.2)  “  {(x^.x^.x^)  en  I  x^+X2=r<R 

6  <  x^  <  1  -  5} 

and  assume  that  R,  5  are  sufficiently  small. 

Let  us  consider  the  edge  AB  in  Figure  3.1.1  as  an  example.  Fig.  3.1.2 
R  S 

shows  the  region  . 


Fig.  3.1.2.  The  edge  neighborhood 
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R  5 

Obviously  the  domain  can  be  described  well  in  the  cylindrical 

coordinates. 

2)  The  neighborhood  of  the  vertex  A  and  the  edge  e  (the  vertex-edge 
neighborhood).  Assume  once  more  that  the  edge  is  given  as  in  (1.1) 

T  ^  ,  2  2  2  2,  ,  2  2  .  , 

Let  X  =  p  (x)  =  x^  +  x^  +  r  (x)  =  +  x^,  sin  <p  =  r/p 

and  denote 

0^'^  =  {x  €  Q  I  p(x)  :s  R,  sin  <p  <  sin 
A 

with  R  and  $  sufficiently  small. 

R  $ 

The  domain  Q  *  associated  to  the  edge  AB  =  e  and  vertex  A  is  shown 
“C  y  A 

in  the  Figure  3.1.3. 

te 


p 

Fig.  3.1.3.  The  vertex-edge  domain 
R  ^ 

Obviously  the  domain  Q  *  can  be  described  well  in  the  spherical 

C-j  A 

coordinates. 

3)  The  neighborhood  of  the  vertex  (not  close  to  the  edges). 
Assume  now  that  in  the  vertex  A  is  the  end  of  the  edges 
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For  example  in  the  Fig.  3.1.1  vertex  A  is  the  end  of  the  edges  =  AB, 

—  _  R 

=  AD,  =  AF.  Denote  now  by  Q  ’  ^  the  edge~vertex  domain  introduced 

R  R  $  1 

above.  We  will  assume  that  $.  are  sufficiently  small  that  f!  ’  i  n  Q  ’  i  = 

1  ^  ej  ,  A  e j ,  A 

0  for  i  *  j.  Then  we  define 

=  {x  €  Q,  p(x)  <  R,  X  ^  i  =  1 . n>. 

A  ,  A 

4)  The  regular  region. 

R  5  R  $  R 

For  the  given  domain  we  defined  the  domains  n  ’  ,  n  ’  ‘  Q  .  Then  we 

•0  ,  A  A 

define  =  £2  -  U  U  -  Un*. 

e  €,A  a 

Let  us  now  make  the  comments  about  selection  of  the  parameters  in  the  4 

particular  domains  we  have  introduced.  We  will  assume  that 

1)  The  domains  of  one  category  (i.e.  edge  or  edge-vertex,  or  vertex 

domains)  do  not  intersect.  On  the  other  hand  we  will  assume  that  the 

parameters  are  selected  to  that  £2^  u  U  £2^‘^‘  u  U  £2^’^!,  uU  £2^  =  £2.  (Note 

°  i  l.A**’*  A  * 

that  the  values  R  are  different  for  different  regions).  The  regions  of 

different  categories  may  intersect  (i.e,  partially  overlap)  so  that  any 

entire  element  J  of  the  usQd  mesh  will  be  in  one  (or  more)  domains  with  one 

exception.  Element  having  the  vertex  coinciding  with  the  vertex  of  the  domain 

Q  —  a  "vertex  element"  has  not  to  be  entirely  in  any  one  region  we 

introduced.  This  exception  is  made  for  practical  reasons  related  to  the  mesh 

R  ( r ) 

generator.  We  see  that  the  domain  £2^  is  analogous  to  the  domain  £2^ 
introduced  in  the  Section  2.1. 

3.2.  The  spaces  and  the  model  problem 

Analogously  as  in  the  Section  2.1  we  introduce  the  spaces  of  functions 
which  norm  is  different  on  every  neighborhood  introduced  in  Section  3.1. 

1)  The  neighborhood  of  the  edge  €  (not  close  to  the  vertex);  the 
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For  example  in  the  Fig.  3.1.1  vertex  A  is  the  end  of  the  edges  =  AB, 

=  AD,  =  AF.  Denote  now  by  the  edge-vertex  domain  introduced 

2  3  “Cj  ,  A 

R  R  $  i 

above.  We  will  assume  that  4>.  are  sufficiently  small  that  Q  ’  *  n  £1  ’  {  = 

1  ,  A  “Cj  ,  A 

0  for  i  ^  J.  Then  we  define 

£2^  =  {x  e  £2,  p(x)  <  R,  x  «  i  =  l,...,n}. 

A  .  , A 

4)  The  regular  region. 

R  5  R  $  R 

For  the  given  domain  we  defined  the  domains  £2  ’  ,  £2  ’  £2. .  Then  we 

e  ,  A  A 


define  £2?  =  £2  -  U  £2^‘^‘-  U  £2^’^1  -  U£2’?. 


ei,A  ^  A' 


Let  us  now  make  the  comments  about  selection  of  the  parameters  in  the  4 

particular  domains  we  have  introduced.  We  will  assume  that 

1)  The  domains  of  one  category  (i.e.  edge  or  edge-vertex,  or  vertex 

domains)  do  not  intersect.  On  the  other  hand  we  will  assume  that  the 

parameters  are  selected  to  that  u  U  u  U  uU  (Note 

0  i  ei  A  ^ 

that  the  values  R  are  different  for  different  regions).  The  regions  of 

different  categories  may  intersect  (i.e.  partially  overlap)  so  that  any 

entire  element  7  of  the  used  mesh  will  be  in  one  (or  more)  domains  with  one 

exception.  Element  having  the  vertex  coinciding  with  the  vertex  of  the  domain 

£2  —  a  "vertex  element"  has  not  to  be  entirely  in  any  one  region  we 

introduced.  This  exception  is  made  for  practical  reasons  related  to  the  mesh 

R  ( r ) 

generator.  We  see  that  the  domain  £2^  is  analogous  to  the  domain  £2^ 


introduced  in  the  Section  2.1. 


3.2.  The  spaces  and  the  model  problem 

Analogously  as  in  the  Section  2.1  we  introduce  the  spaces  of  functions 
which  norm  is  different  on  every  neighborhood  introduced  in  Section  3.1. 

1)  The  neighborhood  of  the  edge  e  (not  close  to  the  vertex):  the 
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R  ^ 

neighborhood  Q  *  .  As  before  we  assume  that  the  edge  e  is  given  in  (1.1) 

R  A  R  5 

and  by  (1.2).  For  given  p,  0  <  p  <  1  we  define  the  space  ) 

R  S 

of  functions  u  on  Q  *  such  that 

e 

— R  5 

i)  u  is  continuous  on  £2^’ 

ii)  for  any  a  =  (a^.a^.a^),  ^  0  integers 

(2.1)  |d“(u(x)  -  u(0,0,X2)|  £  C  d“|r(x) 


(2.2) 


—  u(0,0.X3)|  s  Cd3  a3! 


where  we  denoted 


d  =  (dj,d2,d3),  >  1,  d^  =  d^^d^^d^^ 


a!  =  a^!  a^!  a3! ,  0!  =  1,  |a|  =  |a^+a3+a3l, 


d“u  = 


a'‘‘'u 

a.“>ax“^ax5= 


and  in  (2.1)  and  (2.2)  the  constant  C  is  independent  of  a. 

2)  The  neighborhood  £2  ’  1.  (the  vertex-edge  neighborhood). 

,  A 

Let  X  e  (Xj,X2,X3),  p^(x)  =  ^2  ^  ^3’  “  ^1  ^  ^2’  ^  ~ 

R  $ 

0  <  6,  <  1/2,  0  <  <  1.  Then  by  £_  ^  ’a^  w®  denote  the  space  of  all 

1  ^  PlfRl 

R  $ 

functions  u  on  Q  \  such  that 

e,  A 

-R  $ 

i)  u  is  continuous  on  Q  \ 

e,  A 

ii)  for  any  a  =  ^  integers 


(2.3) 


D  (u(x)  ~  u(0,0,x^) I  ^ 
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(2.4) 


|-^(u(0.0,x.,)  -  u(0, 0,0)1  s  C  p 


•  (j3i+a3-l/2)  ,a3 


where  C  is  independent  of  a. 

R 

3)  The  neighborhood  of  the  vertex  A  (not  close  to  the  edge), 

R 

Let  0  <  p  <  1/2  then  by  £^(£2^)  we  denote  the  space  of  all  functions  u 
R 

on  Q.  such  that 
A 

_R 

i)  u  is  continuous  on  Q, 


ii)  for  any  a  =  (a^.a^.a^),  0  s  integers 


(2.5) 


|d“(u(x)  -  u(0.0,0)|  s  Cp 


and  C  is  independent  of  a. 


4)  The  regular  region  £2^. 

R  R 

By  £o(£2„)  we  denote  the  space  of  all  functions  u  on  £2„  such  that 
p  0  '  u 

— R 

i)  u  is  continuous  on 

ii)  for  any  a  =  (a^,a2,a^),  0  ^  integers 


(2.6) 


|D%!  ^  C  d^a! 


and  C  is  independent  of  a. 

As  was  said  above  we  assume  that  the  regions  are  over taping.  We  define 

the  space  £^(Q)  as  the  space  of  all  u  defined  on  Q  such  that  after 
P 

R  ^  R  ^  R 

constraining  them  on  the  regions  Q  *  ,  Q  '  ,  Q.  and  they  will  belong 

e  e,  A  A  u 

to  £o(£2^’^),  £^  „  £o(f2*)  and  £(£2^).  (We  used  R  which  are 

/3  €  Pi.fe  «.A  P  A  0 

different  in  mentioned  neighborhoods. ) 

Analogously  as  in  Section  1  we  will  consider  the  problem 


(2.7a) 


-Au  =  f  on  £2, 


(2.7b) 


u  =  0  on  r, 


(2.7c) 


au  N„ 

aH  =  8 
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We  will  understand  the  problem  (2,7)  in  the  weak  sense.  Find 

ID  1 

{u  €  H  (Q),  u  =  0  on  r}  such  that  for  any  v  €  H^CQ) 

3 


B(u,v) 


Zdu  dv 

0Xi  dxi 


U  €  hJ(Q)  = 


holds. 

We  will  assume  that  on  every  face  of  the  boundary  of  ^  either  u  =  0 

or  ^  =  g  and  g  is  analytic  on  the  (closed)  face  where  it  is  defined, 
on 

Further  we  will  assume  that  f  is  analytic  on  Q,  In  this  case  we  will  speak 
about  the  problem  (2.7)  with  analytic  input  data.  Then  we  have  proven  in  [18] 


Theorem  2.1,  Let  u  be  the  solution  of  the  problem  (2.7)  with  analytic  input 

data.  Then  u  e  S^(Q)  for  properly  selected  p. 
p 

We  see  that  in  3  dimensions  the  problem  and  the  description  of  the 
regularity  of  the  solution  is  essentially  similar  as  in  two  dimensions  but 
much  more  complex, 

3.3  The  elements 

We  will  assume  that  the  domain  Q  is  covered  by  the  mesh  3JI  which 
partitions  the  domain  Q  into  element  J  and  JR  =  where  J  is  an 

element.  We  will  now  define  elements  and  their  properties  in  various  regions 
we  have  defined  in  Section  3.2.  For  simplicity,  we  consider  only  elements  of 
curved  "brick"  type  i.e.  elements  which  are  images  of  a  standard  cube.  In  a 
similar  way  as  in  Section  2  the  results  hold  for  tetrahedral  which  are  the 

3 

images  of  a  part  of  the  standard  cube.  Let  D  =  I  ,  I  =  (-1,1)  be  the 
standard  cube  with  local  coordinates  17  =  ,  i  =  1,2,3. 

Now  we  describe  the  elements  in  the  various  regions  introduced  in 
Section  2.3. 
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1)  The  elements  in  the  edge  neighborhood  £2 


.R.5 


We  will  consider  two  kinds  of  elements,  ^  and  separately. 

_  R  A 

The  element  J.  Let  ?  c  £2  ’  (i.e.  distance  of  3"  to  the  edge  e  is 

-  e 

_  3 

positive)  be  the  image  of  D,  D  =  I  ,  I  =  (-1 , 1)  1  by  a  mapping  Mg..  Assume 


that  Mg.  is  a  one  to  one  mapping  of  D  onto  S'  with 

S'  =  {x  I  =  XAn),  7>  e  D,  i  =  1,2,3} 


and  X^Ct))  are  analytic  functions  on  D. 


Denote 

(3.1a) 


kO')  =  min  r  (x) , 
xeJ 


(3.1b) 


K(y)  =  majc  r  (x) . 
xeS 


We  assume  about  ff  and  M  the  following: 
i)  There  exists  constant  1  <  A  <  oo  such  that 


(3.2a) 

k(3) 

ic(3) 

:£  A 

ii) 

For 

|a|  =  m,  m  >  0,  m 

an  integer  and  any  7)  €  D 

(3.2b) 

|d“x.(7,)| 

B  O 

O 

U 

VI 

m!  K,  i  =  1 , 2, 

(3.2c) 

\D^x^iin)\ 

<  C  d"* 
“  0  0 

m!  H,  H  <  Hq 

with  Cq 

and 

H  independent  of 

m  and 

K  (and 

^  such  that 

iii) 

The 

Jacobian  determinant 

|J| 

satisfies 

Cj  H  s 

|J|  = 

a(Xi,X2,X3) 
9(t)i.D2  Vs) 

(3. 2d) 

where  and  are  independent  of  k  and  H. 


37 


Obviously  A  >  1  in  (3.2a).  Using  (3.2a-d)  we  easily  see  that  also 


(3.2e) 


kO") 

k(J) 


a  A* 


where  A*  >  1  depending  on  constants  in  (3.2b,c,d).  Let  us  note  that  the 
constants  A,  A*,  C^,  d^,  C^,  C^,  are  not  mutually  independent;  there  are 
relations  between  them.  Further  we  assume  that  these  constants  are  the  same 
for  all  elements  of  the  meshes  3H  e 


2)  The  element  ^  .  The  element  J  is  the  image  of  D  by  the  mapping 

-  e  ^  J 

satisfying  the  conditions  (3.2b,c,d)  and  with  one  edge  of  3^^  on  the  axis 
which  is  a  part  of  the 'edge  e  of  the  domain  Q  i.e.  F  n  9'^  ^  0.  For 
concreteness  and  without  any  loss  of  generality  we  assume  that 

-1,7)3)  =  X^C-l, -1,7)3)  "  ^ 

Remark  3.1;  We  see  from  (3.2a)(3.2e)  that  the  size  of  the  element  J  in 
(x^,X2)  is  of  the  order  of  the  distance  of  J  from  the  origin  i.e.  the  edge 
-e  and  the  size  of  the  element  in  the  variable  X3  is  of  order  H.  This 
means  that  the  element  J  is  a  "needle"  i.e.  has  a  very  large  aspect  ratio. 

In  the  variables  ^^>^2  element  S'  has  the  same  character  as  in  the  two 

dimensional  case. 

R  ^ 

2)  The  elements  in  the  vertex-edge  neighborhood  Q  ’  . 

^  y  ri 

1)  The  element  J.  Let  W  c  (i.e.  W  r\  T  =  0)  be  the  image  of  D 

-  e,  A 

by  an  analytic  mapping  M^..  We  assume  that 

=  Xj(7)),  7)  e  D,  i  =  1,2,3} 

where  XAti)  are  analytic  functions  with  the  properties  spelled  out  below. 
Let 


I 
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(3.4a) 

K. (J)  =  min  p(x) , 
xeJ 

(3.4b) 

K,  O')  =  max  p(x) , 

xeS" 

(3.4c) 

K  (?)  =  min  sin  <p(x), 

X€f 

(3.4d) 

K„(?)  =  max  sin  ^(x). 
^  xe? 

About  the  element  J  and  the  mapping  we  will  assume  the  following 

i)  There  exist  constants  1  <  A^  <  «  such  that 


(3.5a) 


am  <  A 

ica(J)  ■ 


(3.5b) 


.  A 

KzO-)  “  2 


ii)  For  |a|  =  m,  m  >  0  an  integer 


(3.5c) 


|d“x^(t))|  s  Cpd™  m!  i  =  1,2 


(3.5d) 


|d“X3(t})1  £  CpdJ  m! 


with  Cq  independent  of  m.ic^^.K^. 

iii)  The  Jacobian  determinant  |J|  satisfies 


(3.5e) 


a(Xi.X2,X3)  ^  ^  3  2 

d{-ni.'n2>'n3'i  i  i  2 


where  C. 

1 

Combining 

(3.5f) 


and  C. 
1 

(3.5a-e) 


are  independent  of 

m 

there  also  exist 


KiCJ) 

Kim 


2' 

so  that 


* 


2:  A. 


(3.5g) 


<2  0-)  > 
icaO')  “ 


» 
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where  A^,  >  1  depend  on  the  constants  in  (3.5a-e). 

Remark  3.2.  We  see  from  (3.5a)  (3.5b)  that  the  element  proportions  depend  on 
Kj  and  and  its  size  is  of  the  order  of  the  distance  to  the  origin  and 

axis  x^. 

The  element  3"  ..  We  assume  that  3  .  is  image  of  D  by  the  mapping 

satisfying  the  properties  (3.5c-e)  but  instead  (3.5ab)  we  will  assume  only 

(3.5a)  and  that  one  edge  of  3"  lies  on  the  axis  x  which  coincides  with 

the  edge  e  i.e.  W  n  T  *  <t>  but  A  3  .  For  concreteness  we  will 

0  j  A  0/\ 

assume  that  (-1 , -1  .t)^)  =  X2(“1.“1.t12^  ” 

Remark  3.3.  We  described  only  the  elements  for  which  A  i  ?. 

Elements  3"  for  which  A  e  J  will  be  addressed  in  the  subsection  5  below. 

3)  The  elements  in  the  region  £^. 

_  p 

Let  3"  c  n.  be  as  before  the  image  of  D  so  that 
A 

y  =  (x  I  Xj^  =  Xj^('r>),  T)  6  D,  i  =  1,2,3} 

and  X^(7})  are  analytic  functions  with  the  following  properties.  Let 

(3.6a)  ic(3')  =  min  p(x), 

X63' 


(3.6b) 


k(3)  =  max  p(x) . 
X€3 


Then  we  will  assume  the  following: 

i)  There  exists  a  constant  1  <  A  <  oo  such  that 

k(3-) 


(3.7a) 


ic(3) 

ii)  For  |a|  =  m,  m  >  0  integer 


£  A. 


(3.7b) 


IdVi  i  Cpd™  m!  K,  i  =  1,2,3 
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with  Cq  independent  of  m  and  ic. 

iii)  The  Jacobian  determinant  |J|  satisfies 


(3.7c) 


Ck  5  |J|  = 


a(Xi,X2,X3) 

d(Vl,V2,V3) 


where  C  and  C  are  independent  of  k.  Combining  (3.7abc)  we  also  get 


(3.7d) 


ic(J) 

kO-) 


* 

^  A 


where  A  >  1  depending  on  the  constants  in  (3.7abc). 


Remark  3.4.  We  see  that  the  element  has  not  large  aspect  ratio  and  its  size 
depend  on  the  distance  from  the  vertex  A. 


Remark  3.5.  Here  we  do  not  address  the  element  J  with  A  e  it  will  be 
addressed  in  the  subsection  5  below. 

4)  The  elements  in  the  regular  region  Qq. 

T5 

Once  more  let  7  e  £2^  and  3"  be  image  of  D  by  the  mapping  =  {X^^} 
i  =  1,2,3.  Let 

(3.8)  K  =  max  (dist  x,r). 

xeJ 

i)  We  will  assume  th^t  for  |al  =  m,  m  >  0  integer 


(3.9a) 


|dV(t?)|  s  C^dJ  m!  K 


with  Cq  independent  of  m  and  k. 

ii)  The  Jacobian  determinant  |J|  satisfies 


(3.9b) 


Ck  s  |J|  = 


a(Xi,X2,X3) 

d(vi,V2.V3^ 


£  Ck" 


where  C  and  C  are  independent  of  k. 


41 


*  R 

From  (3.9a,b)  it  follows  that  the  volume  of  any  U  €  is  bounded  from  below 
and  so  only  finite  number  of  elements  S'  e  exist  and  k  is  equivalent  to 
the  diameter  of  J. 

R 

Remark  3.6.  All  elements  in  have  no  large  aspect  ratio. 

5)  The  vertex  element  S'.. 

A 

So  far  we  did  not  address  the  element  S"^  such  A  e  In  all  previous 

cases  we  have  assumed  that  the  (entire)  element  S'  is  in  one  (or  more) 
domains.  About  the  vertex  element  we  will  assume  that  in  general  ST  is  not 

R  ^  R  - 

necessarily  in  the  interior  of  any  of  the  domains  A  €  S'.  We 

will  assume  that  S*^  is  the  image  of  D  by  the  mapping  =  {X^}i  =  1,2,3. 
Denoting  diam  S'  =  k  then ’we  will  assume  that  (3.5  bcde)  (3.7  ab)  hold. 
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3.4  The  mesh  and  the  finite  elements 


We  will  consider  a  family  3^  of  meshes  on  n.  A  mesh  OH  €  ?  will  be  a 

partition  of  £2  into  the  set  of  elements  J  e  M,  U  =  £2.  We  will  assume 

yeOJl 

that  the  elements  3"  are  curvilinear  bricks  which  properties  were  described 
in  the  Section  3.3.  For  simplicity  we  are  restricting  ourselves  to  the  case 
of  brick  elements  only,  although  the  theory  holds  for  tetrahedrons  in  the 
similar  way  as  in  two  dimensions  the  triangular  elements  were  treated. 

To  every  ?  €  OH  we  associate  an  analytic  mapping  which  is  defined 

on  the  standard  cube  and  7  is  the  image  of  the  cube  by  the  mapping  M^..  In 
an  obvious  way  we  can  speak  about  the  vertices,  edges  or  faces  of  the  elements 


S'  €  511. 

As  usually  we  will  assume 

that 

if 

^i’ 

?.  e  DH 

J 

then 

either 

y .  n  S' .  = 
1  J 

(f>  or 

7.  n  7  .  is  a 

1  J 

common  vertex 

of 

J. 

1 

and 

S’ .  or 

J 

W .  is 
J 

a  (common) 

entire 

edge  or  entire 

face  of 

1 

and 

S'.. 

J 

If 

S’. ,3-. 

1  J 

have 

common 

edge  or 

face  we  assume  that  the  mapping  and  have  the  usual  properties 

Ji  J  j 

guaranteeing  the  continuity  of  the  finite  elements.  The  meshes  consisting  of 
elements  described  above  have  a  special  character  of  a  geometrical  mesh  in  the 
neighborhood  of  edges  and  vertices.  This  follows  from  (3.2a)  (3.2e),  (3. Sab), 
(3.5fg)  (3.7a)  (3.7d). 

We  will  consider  the  family  3  of  meshes  characterized  by  two  parameters 

(<r,n)  analogously  as  in  two  dimensions.  Consider  the  edge  e  =  {x^  =  =  0, 

0  <  x^  <  1}  and  the  vertex  A  =  (0,0,0). 

R  5 

1)  The  mesh  in  the  edge  region  Q^’  .  Here  the  mesh  is  geometric  in 
(Xj,X2)  with 

Kitf  )  ^  Co''^. 
e 

In  the  Xg  variable  the  series  of  the  elements  are  independent  of  n, 
i.e.,  the  mesh  is  fixed  in  the  x^  direction;  see  (3.1),  (3.2).  We  can 
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(analogously  as  in  two  dimensional  case)  speak  about  the  layer  (in 

and  the  level  which  is  a  sequential  number  of  element  in  direction.  The 

mesh  hence  has  n-layers  and  independent  of  n)  levels.  Any  layer  has  at 

R  6 

most  K  elements  and  hence  the  total  number  of  elements  in  Q  '  is  0(n) 

e 

similarly  as  in  2  dimensions. 

R  $ 

2)  The  mesh  in  the  vertex-edge  region  ^ 

The  mesh  has  geometric  character  in  (x^.x^)  and  also  in  x^.  We  can 
also  speak  here  about  the  layers  (in  and  levels  in  x^.  The  mesh  (see 

(3.4),  (3.5))  is  such  that 


,)  £  Co'” 

2  -e,  A 


and 


min 


R,$ 
e.y  A 


K^O")  ^  Qr^. 


By  the  analogous  argument  as  in  the  two  dimenional  settings  we  can  see  that 

the  number  of  elements  in  is  of  order  O(n^). 

e,  A 

R 

3)  The  mesh  in  the  vertex  region 

Here  the  mesh  is  geometric  in  the  p(x)  variable  and  quasiuniform  in  the 
two  other  variables  (see  (3.6),  (3.7))  with 


Co-”  £  min  kO")  £  Qr^. 

A 

R  $ 

Hence  the  number  of  elements  in  Q  * .  is  of  order  0(n). 

e,  A 

p 

4)  The  mesh  in  the  regular  region  £1^. 

Here  the  mesh  has  number  of  elements  independent  of  n.  Hence  the  number 
of  elements  in  is  of  order  0(1). 

5)  The  set  of  vertex  elements. 

Here  we  have 
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Co''^  £  diam  O^q)  -  Co'^. 

Obviously  we  have  0(1)  vertex  elements. 

Remark  4.1.  Similarly  as  in  2  dimensions  (see  remarks  4.1  and  4.2  in  the 

Section  2)  we  have  analogous  inequalities  for  k(3^).  For  example  let  the 

element  7  c  be  located  in  the  layer  and  level.  Then 

e,  A 

C<r  pj  s  k^(7)  s  Co-  P2  . 

and 

Co-”p"^  s  »c^(J)  s  C0%"^ 

with  0  <  Pj.p^  <  1. 

In  practice  we  construct  the  mesh  with  ~  ~  number  of  elements 

R  $ 

in  the  region  Q  \  can  be  estimated  analogously  as  in  the  two  dimensional 
case.  Similar  estimates  for  the  elements  located  in  other  regions  are 
analogous . 

The  shape  functions  are  given  on  the  master  elements  which  in  our  case  is 
the  unit  cube  D.  We  will  assume  that  the  elements  are  of  degree  pO")  in 
every  variable  i  =  1,2,3  separately.  For  simplicity  we  will  assume  that 

pO")  =  p  i.e.  the  degrees  are  uniform  although  we  could  consider  pO") 
different  for  every  element  and  different  in  every  direction.  This  would 
increase  the  effectiveness  of  the  method.  The  space  of  the  finite  elements  is 
denoted  by  S(3J),p)  =  5(0-, n,p)  c  Hq(£2). 

3.5.  The  rate  of  convergence  of  the  h-p  version. 

In  the  Section  3.4  we  have  described  the  family  of  the  meshes  7 
characterized  by  two  parameters  ((r,n).  Given  the  exact  solution  u^  with 
the  properties  given  in  the  section  3.2  we  follow  the  framework  outlined  in 
the  Section  2.7.  First  we  construct  the  approximation  element  by  element 
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and  then  the  corrections  which  guarantee  the  conformity  of  the  elements. 

By  a  detailed  analysis  of  the  approximation  of  the  functions  u  €  £(£2) 
we  can  show  (see  [16])  that  for 


pO")  =  p  =  fxn,  fi  >  0 

R  5  R  $  R 

there  exists  a  function  x  ^  S(<r,n,p)  so  that  in  any  regions 


Qq  we  have 
(5.1) 


where  r  >  0  depends  on  the  constants  in  introduced  in 

Section  3.2. 

Let  us  now  estimate  the  error  in  the  term  of  degrees  of  freedom. 

1)  The  region 

R  3 

As  we  have  seen  above  there  are  0(n)  elements  in  the  region  S2^’  . 

o  R  5 

Every  element  has  0(p  )  degrees  of  freedom  NO”).  Denoting  by  )  the 

R  S 

number  of  degrees  of  freedom  in  we  get 

N(£2^’^)  =  O(np^)  =  0(n'^) 
when  we  used  pCf)  =  pn.  Hence  we  have 


(5.1) 


I  ^  I  .1  R  ^ 

e 


Or 


-m1/4 

rN 


=  Ce 


1/4 


2)  The  region 

e  >  A. 

R  $  2 

As  stated  above  the  total  number  of  elements  in  ^  \  is  0(n  ). 

e  j  A 


Hence 


N(n^’t)  =  O(n^p^)  =  O(n^)  and 
■e,  A 


(5.2) 


,  ,  ^  „  -r*N 

e,  A 


1/5 


3) 


The  region  £2^. 

In  n.  there  are  0(n)  elements  and  hence 
A 
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(5.3) 


-1/4  *.,1/4 

I  V^l  1  R  "  ^  ■ 


4)  The  region 

The  number  of  elements  in  this  region  is  0(1)  and  hence 

-.,1/3  *„l/3 

(5.4)  1  R  “  =  Ce  ^ 

°  eU^I) 

5)  The  vertex  elements. 


Ce 


Here  we  have  only  0(1)  elements  and  hence  the  error  is  of  order 


Hence  we  have  proven 


Theorem  3.5.1.  Let  (see  Section  3.2)  be  the  exact  solution  of 

the  problem  (2.7).  Further  let  the  meshes  JIKcTjO),  <r  <  1  are  as  in  the 
Section  3.4  and  let  p(J)  =  pn,  ^  >  0,  properly  chosen.  Then 

(5.5)  llu„  -  u,,,,!,,  s  C  inf  1“  -  >:l  l  “ 

;^€S(cr,if),p)  H  (Q) 

where  N  =  dim  S(<r,7j,p)  is  the  number  of  degrees  of  freedom,  y  and  C 

depend  on  £^(Q),  the  distortion  of  the  elements  the  solution  u^  and  the 
p  u 

domain  Q  but  are  independent  of  N. 


Remark  5.1.  Theorem  3.5.1  follows  from  (5. 2) (5. 3) (5. 4) .  We  see  that  the 
1/5 

factor  N  is  due  to  the  vertex-edge  singularly  of  the  solution.  Hence  we 

can  expect  that  the  rate  e  will  be  visible  only  for  large  N  i.e. 

-1 

— vfj  1 

high  accuracy  and  for  smaller  accuracies  we  can  expect  the  rate  e  ,  ^<  a  < 
1 

We  will  see  it  in  the  next  section. 

Remark  5.2.  We  assumed  in  the  Theorem  5.1  that  the  (pull-back)  polynomials 
are  at  the  same  degree  p  in  elements.  It  is  possible  to 
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prove  that  the  error  in  (5.5)  holds  with  better  constants  if  the  degree  of  the 
elements  is  different  in  different  elements  and  in  different  directions  in 
i  =  1,2,3. 

Remark  5.3.  In  [12]  we  have  proven  that  in  the  one  dimension  and  the 

function  of  the  analogous  type  as  here  the  exponential  rate  is  c  .  a  =  1/2 

and  a  cannot  be  made  larger  for  any  mesh  and  any  degrees  of  the  polynomials. 

“"TN  1 

We  conjecture  that  the  exponential  rate  e  ,  a  =  ^  in  3  dimensions  cannot 
be  improved  also  if  any  mesh  and  any  anisotropic  degree  distributions  would  be 
considered. 

Similarly  we  can  prove 

Theorem  3.5.2.  Let  u^  be  the  exact  solution  of  the  Problem  2.7  and  be 
analytic  on  Q,  Consider  the  meshes  !IH((r,7})  with  <r  =  1  (hence  3Il((r,7j)  has 
only  finite  number  of  elements  independently  of  n)  and  let  pO”)  =  pn,  p  >  0. 
Then 

y  1/3 

||u  -  U  II  i  C  inf  |u  -  xl  ^  Ce"^”  . 

X6S((r,7j,p) 

Remark  5.4.  For  the  results  mentioned  above  see  also  [17]. 

3.6.  A  numerical  example 

In  this  section  we  will  show  a  typical  numerical  example.  Consider 
domain  shown  in  the  Fig.  3,6.1  and  3.6.2  and  the  problem 

Au  =  0  on  Q. 

The  boundary  conditions  are 

i)  u  =  0  on  BFHGQKD,  ABDC,  AEFB,  ACNE,  PRKQ,  SGHF,  BFHGQKD, 

ii)  =  0  on  SLRPRQG,  MJHFEN,  MNCDKRL, 

on  the  bottom  surface. 
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Fig.  3.6.1 


The  domain  n  and  the  scheme  of  the  mesh 


We  constructed  the  geometrical  meshes  of  the  type  described  in  the  previous 
subsections  with  <r  =  0.15.  One  of  the  scheme  of  the  meshes  (elements 


surfaces)  is  shown  in  the  figures.  Analogously  as  in  the  two  dimensional  case 
for  every  mesh  characterized  by  n=i,  i=l,...,7  we  solve  the  problem  by 
p  =  2 . 8.  Theorem  3.5.1  suggests  that  the  error  measured  in  the  energy 

norms  behaves  as  e  "  where  a  s  5  and  we  can  expect  that  we  will  see 
4  <  a  £  5  in  the  computation.  Hence  we  plot  the  relative  error 
the  scale  log  llelli_  x  Then  we  will  expect  that  of  the  error  of 

the  h-p  version  will  decay  linearly.  In  the  Fig.  3.7.3  we  show  the  errors  for 
a  =  4  and  in  Fig.  3.7.4  for  a  =  5. 


Fig.  3.6.3  The  error  as  function  of  the  degrees  of  freedom  in  the 

scale  log  x 
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Fig.  3.6.3  The  error  as  function  of  the  degrees  of  freedom  in  the 

scale  log  llelljj^  x 

We  see  that  the  h-p  version  converges  exponentially  with  ** 

Ce  4  <  a  ^  5  as  expected  in  the  range  of  the  computation. 

Remark  6.1.  The  exact  solution  is  of  course  not  known.  We  computed  the 

strain  energy  &  of  the  exact  solution  by  extrapolation  and  the  error  of  the 

2 

finite  solution  is  then  ||e||„  =  ^(uexact)  ”  ^(ufe)- 

jb 

Remark  6.2.  The  above  example  was  computed  by  Dr.  B.  Anderson  (Aeronautical 
Institute  of  Sweden)  using  the  program  STRIPE  (developed  in  the  Aeronautical 
Institute) . 


51 


References 


1.  Babuska,  I.,  Guo,  B.Q.  (1987).  The  Theory  and  Practice  of  the  h-p 

version  of  the  Finite  Element  Method.  Comp.  Meth.  in  Partial  Diff. 
Equations  -  VI,  R.  Vichnevesky  R.S.  Stepleman  eds. ,  241-247. 

2.  Babuska,  I.,  Suri,  M.  (1994).  The  p  and  h-p  versions  of  the  finite 
element  method.  An  overview.  Comp.  Meth.  Appl.  Mech.  Engg.  80,  5-26. 

3.  Babuska,  I.,  Guo,  B.Q.  (1992).  The  h,p  and  h-p  version  of  the  finite 

element  method;  basic  theory  and  applications.  Advances  in  Engineering 
Software,  159-174. 

4.  Babuska,  I.,  Guo,  B.Q.  (1988).  Regularity  of  the  solutions  of  elliptic 

problems  with  piecewise  analytic  data.  Part  I.  Boxmdary  Value  Problems 
for  Linear  Elliptic  Equations  of  Second  Order,  SIAM  J.  Math.  Anal.  19, 
172-203. 

5.  Babuska,  I.,  Guo,  B.Q.  (1989).  Regularity  of  the  solutions  of  elliptic 

problems  with  piecewise  analytic  data.  Part  II.  The  trace  spaces  and 
applications  to  the  Boundary  value  problem  with  Non-homogeneous  boundary 
conditions.  SIAM  J.  Math.  Anal.  20,  765-781. 

6.  Guo,  B.Q.,  Babuska,  I.  (1993).  On  the  regularity  of  elasticity  problems 

with  piecewise  analytic  data.  Advances  in  Applied  Mathematics  14, 
307-347. 

7.  Babuska,  I.,  Guo,  B.Q.,  and  Osborn,  J.E.  (1989).  Regularity  and 

numerical  solutions  of  eigenvalue  problems  with  peicewise  analytic  data, 

SIAM  J.  Num.  Anal.  26,  1534-1560. 

8.  Babuska,  I.,  Guo,  B.Q.  (1989).  The  h-p  version  of  the  finite  element 

method  for  problem  with  non-homogenous  essential  boundary  conditions. 
Comp.  Meth.  Appl.  Math-  Engg.  74,  1-28. 

9.  Guo,  B.Q.,  Babuska,  I  (1986).  The  h-p  version  of  the  finite  element 

method.  Part  I.  The  basic  approximation  results.  Comp.  Mech.  1,  21-41. 

10.  Guo,  B.Q.,  Babuska,  I.  (1986).  The  h-p  version  of  the  finite  element 
method.  Part  II.  General  results  and  applications.  Comp.  Mech.  1, 
203-220. 

11.  Babuska,  I.,  Guo,  B.Q.  (1988).  The  h-p  version  of  the  finite  element 
method  for  domains  with  curved  boundaries.  SIAM  J.  Numer.  Anal.  25, 
837-861 . 

12.  Gui,  A.,  Babuska,  I.  (1986).  The  h-p  versions  of  the  finite  element 
method  in  one  dimension.  Part  I,  the  error  analysis  of  the  p-version; 
Part  II,  the  error  analysis  of  the  h  and  h-p  versions,  Numer.  Math.  43, 
577-612,  613-657. 


52 


13.  Guo,  B.Q.  (1988).  The  h-p  version  of  the  finite  element  method  for 
elliptic  equations  of  2m  order.  Num.  Math.  53,  199-224. 

14.  Babuska,  I,  Eleman.  H.C. ,  Markley,  K.  (1992).  Parallel  implementation  of 
the  h-p  version  of  the  finite  element  method  for  a  shared  memory 
architecture  SIAM  J.  Sci.  STat.  Comp.  13,  1433-1450. 

15.  Babuska,  I.,  Elman,  H.C.  (1993).  Performance  of  the  h-p  version  of  the 
finite  element  method  with  various  elements.  Internat.  J.  Num.  Meth. 
Engg.  36,  2503-2533. 

16.  Babuska,  I.,  Guo,  B.Q.  The  h-p  version  of  the  fintie  element  method  in  3 
dimenions.  To  appear. 

17.  Guo,  B.Q.  (1994).  The  h-p  version  of  the  finite  element  method  for 
solving  boundary  value  problems  in  polyhedral  domains.  Boundary  value 
problems  and  integral  equations  in  nonsmooth  domains.  M.  Costabel, 

M.  Dauge,  S.  Nicaise,  eds. ,  101-120. 


53 


The  Laboratory  for  Numerical  Analysis  is  an  integral  part  of  the  Institute  for  Physical 
Science  and  Technology  of  the  University  of  Maryland,  under  the  general  administration  of  the 
Director,  Institute  for  Physical  Science  and  Technology.  It  has  the  following  goals: 

To  conduct  research  in  the  mathematical  theory  and  computational  implementation  of 
numerical  analysis  and  related  topics,  with  emphasis  on  the  numerical  treatment  of 
linear  and  nonlinear  differential  equations  and  problems  in  linear  and  nonlinear  algebra. 

To  help  bridge  gaps  between  computational  directions  in  engineering,  physics,  etc.,  and 
those  in  the  mathematical  community. 

To  provide  a  limited  consulting  service  in  all  areas  of  numerical  mathematics  to  the 
University  as  a  whole,  and  also  to  government  agencies  and  industries  in  the  State  of 
Maryland  and  the  Washington  Metropolitan  area. 

To  assist  with  the  education  of  numerical  analysts,  especially  at  the  postdoctoral  level, 
in  conjunction  with  the  Interdisciplinary  Applied  Mathematics  Program  and  the 
programs  of  the  Mathematics  and  Computer  Science  Departments.  This  includes  active 
collaboration  with  government  agencies  such  as  the  National  Institute  of  Standards  and 
Technology. 

To  be  an  international  center  of  study  and  research  for  foreign  students  in  numerical 
mathematics  who  are  supported  by  foreign  governments  or  exchange  agencies 
(Fulbright,  etc.). 

Further  information  may  be  obtained  from  Professor  I.  BabuSka,  Chairman,  Laboratory  for 
Numerical  Analysis,  Institute  for  Physical  Science  and  Technology,  University  of  Maryland,  College 
Pmk,  Maryland  20742-2431. 


